Quantum dynamics of Bose-Hubbard Hamiltonians beyond Hartree-Fock-Bogoliubov: 

The Bogoliubov backreaction approximation 
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We formulate a method for studying the quantum field dynamics of ultracold Bose gases con- 
fined within optical lattice potentials, within the lowest Bloch-band Bose-Hubbard model. Our 
formalism extends the two-sites results of Phys. Rev. Lett. 86, 000568 (2001) to the general 
case of M lattice sites. The methodology is based on mapping the Bose-Hubbard Hamiltonian to an 
SU(M) pseudospin problem and truncating the resulting hierarchy of dynamical equations for corre- 
lation functions, up to pair-correlations between SU(M) generators. Agreement with few-site exact 
many-particle calculations is consistently better than the corresponding Hartree-Fock-Bogoliubov 
approximation. Moreover, our approximation compares favorably with a more elaborate two-particle 
irreducible effective action formalism, at a fraction of the analytic and numerical effort. 

PACS numbers: 3.75.Kk 



I. INTRODUCTION 

Strong correlation effects, which imply enhanced quan- 
tum fluctuations around mean field order parameters, are 
playing an increasingly important role in recent experi- 
ments on dilute quantum gases. One strategy for boost- 
ing the importance of correlations and fluctuations in- 
volves the control of coupling parameters. Interatomic 
interactions can be effectively tuned by means of mag- 
netic Feshbach resonances 0, 0, 0, 13 IE 0- allowing 
for a controlled transition into the non-unitary regime 
n 1//3 a s > 1, where the effective s-wave scattering length 
a s is larger than the average distance between particles 
n 1 / 3 with n being the number-density of the gas. Quan- 
tum fluctuations also dominate quasi-one-dimensional 
systems where trans- 

verse confinement may be used to increase the effective 
coupling strength g\jj = 2h 2 a s / '(mlj_) without explicit 
control of the three-dimensional s-wave scattering length. 
In the extreme Tonks- Girardeau strong-coupling regime 
gi£>m/ (h 2 n) 3> 1, spatial correlations dictate the impen- 
trability of bosons, leading to ideal fermion like density 
distributions 

An alternative to increasing effective interaction 
strengths, is to decrease other (e.g. kinetic) terms in 
the many-body Hamiltonian. In a Bose gas confined by 
an optical lattice, an effective momentum cutoff is intro- 
duced by controling the barrier heights, thus suppressing 
the hopping frequency J between adjacent sites. Given 
N particles interacting with strength U, the strong- 
interaction regime is achieved for UN/ J > 1, as man- 
ifested in the quantum transition from a superfluid to a 
Mott-insulator phase (lj ITl ITs| . 

Considerable theoretical effort is currently aimed at 
developing efficient methods for the description of cor- 
related quantum gases far from equilibrium. One ap- 
proach relies on perturbations of the lowest-order mean- 
field theory given by the Gross-Pitaevskii (GP) equation. 
The result is a family of mean-field pairing theories. The 



standard zero-temperature Bogoliubov prescription |l9j 
gives the natural small-oscillation modes by lineariza- 
tion about the GP ground-state. However, this linear 
response theory does not account for the backreaction 
of excitations on the condensate order-parameter and is 
thus limited to small perturbations and short timescales. 
Backreaction is accounted for within the Hartree-Fock- 
Bogoliubov (HFB) theory, which prescribes a set of cou- 
pled equations for the condensate order-parameter and 
pair correlation functions [2pl l2lt l22l |23| . Since both 
normal and anomalous correlations are included, this ap- 
proach comes at the cost of ultraviolet divergences of 
anomalous quantities. While this problem is relatively 
easy to deal with by renormalization of the coupling pa- 
rameters, a more serious issue, also related to the inclu- 
sion of anomalous correlations, is the HFB spectral gap 
|2dj |. This unphysical gap in the excitation spectrum re- 
sults in from the breaking of U{\) gauge symmetry and 
the consequent elimination of the Goldstone modes cor- 
responding to gauge transformations of the broken sym- 
metry solution. An intermediate form between Bogoli- 
ubov and HFB is the HFB-Popov (HFB-P) approxima- 
tion |2Ctl24| where U(l) symmetry is restored by elimina- 
tion of noncondensate anomalous terms only. While the 
resulting theory is gapless, it does not conserve the total 
number of particles and is thus inadequate for describing 
dynamical condensate depletion. Finally, if all anomalous 
quantities are neglected, one obtains the bosonic Hartree 
Fock (HF) theory [2(j which is both gapless and conserv- 
ing, but does not allow for any dynamical depletion, since 
the populations of condensed and non-condensed parti- 
cles are conserved separately. It is thus highly desirable 
to develop a theoretical description that (a) is U(l) in- 
variant and hence gapless, (b) conserves the total number 
of particles, yet (c) allows for dynamical depletion of the 
condensate. 

Recently, a perturbative approximation scheme based 
on a two-particle irreducible (2PI) effective action ex- 
pansion, has been used to study the noneq uilibrium dy- 
namics of condenstaes in optical lattices [23| within the 
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lowest Bloch-band Bose-Hubbard model. Within the 2PI 
effective action expansion, the Bogoliubov and HFB the- 
ories emerge as one-loop and a single two-loop correc- 
tion respectively, to the classical GP action. Higher-order 
approximations, obtained by including two-vertex terms 
in the diagramatic expansion of the effective action (de- 
noted as in Ref. (2^1 by '2nd') and by a 1 /N expansion up 
to second-order in the coupling strength (denoted hence- 
forth by T/7V') with Af being the number of auxiliary 
classical fields used to approximate the quantum-field, 
have been compared with HFB and exact few-sites nu- 
merical calculations. The results demonstrate some im- 
provement of the higher-order approximations over HFB 
in predicting the exact many-body dynamics. However, 
at sufficiently long times all approximations fail due to 
interaction effects. A nonperturbative 1/Af 2PI effective 
action expansion approach have also been developed and 
applied to the equilibration of a homogeneous Bose gas 
in ID El. 



In this work we develop a mean-field theory for the 
description of quantum dynamics in the Bose-Hubbard 
model. The technique, referred to here as Bogoliubov 
Back Reaction (BBR), is a many-site extension of previ- 
ous work on a two-site model |2y, , based on the per- 
turbation of equations of motion for the reduced single- 
particle density operator, instead of the usual field oper- 
ator approach. The resulting equations involve the two- 
point reduced single-particle density matrix (SPDM) and 
the four-point correlation functions. They contain only 
normal (i.e. number conserving) quantities, and are thus 
U{1) symmetric. The approximation conserves the total 
number of particles, yet it allows for population trans- 
fer from the condensate to the excitations, thus account- 
ing for condensate depletion during the evolution. We 
compare BBR calculations with full many-body numer- 
ical results for up to a hundred particles and five lat- 
tice sites, as well as with HFB and 2PI effective action 
results. The BBR results give better, longer-time pre- 
dictions than current rival approximations, at a small 
fraction of the theoretical effort. 



In section II we present the Bose-Hubbard model and 
the standard HFB approach. In section III we trans- 
form the Bose-Hubbard Hamiltonian with M lattice sites 
into an SU (M) pseudospin problem, derive dynamical 
equations for the SU (M) generators spanning the single- 
particle density operator, and truncate the resulting hier- 
archy of dynamical equations for correlation functions to 
obtain the BBR equations of motion. Section IV contains 
numerical few-sites results and comparison with HFB as 
well as 2PI effective action approximation methods. Dis- 
cussion, conclusions and prospects for future research are 
presented in section V. 
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FIG. 1: (color online) Population imbalance w in a two-site 
system as a function of rescaled time r for L = 2 (a), L = 
4 (b), L = 5 (c), and L = 10 (d). The total number of 
particles is set to N = 100. Solid blue lines, corresponding to 
exact many-body numerical results, are compared to the GP 
(dotted lines), HFB (dash-dotted lines), and BBR (dashed 
lines) approximations. 



II. CONVENTIONAL MEAN-FIELD THEORIES: 
GROSS-PITAEVSKII AND 
HARTREE-FOCK-BOGOLIUBOV 



We begin with the standard Bose-Hubbard model 
Hamiltonian for an ultracold gas in a one-dimensional 
periodic optical lattice 

H = (4fiOi + aja i+ i) a K T a^ (1) 



where di and aj are annihilation and creation operators 
respectively, for a particle in site i. We consider only on- 
site interactions with strength U and nearest-neighbor 
tunneling with hopping rate J. These approximations 
are justified because adjacent site interactions and next- 
to-nearest-neigbhbor tunneling amplitudes are character- 
istically at least two orders of magnitude smaller than 
on-site interactions and nearest-neighbor hopping |16| . 
The Bose-Hubbard model is viable as long as there are 
no transitions into excited Bloch bands. 

Using the Hamiltonian we write the Heisenberg 
equations of motion for the field-operators dj , 



a i 
i— dj — J [a,j-x + &j+i) + Uajdjdj 



(2) 



The lowest order mean field theory for the Bose-Hubbard 
model is obtained by replacing the field operators dj , and 
dj the by c- numbers aj and a*. This approximation is 
tantamount to assuming coherent many-body states with 
a well-defined phase between sites. Rescaling a — > y/Na 
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and r = Jt we arrive at the discrete GP equation, 

l ^. a i = ( a i-i + a J+i) + L M 2a i > ( 3 ) 

where L = UN/ J is the characteristic coupling parame- 
ter. Within the GP mean field theory © fluctuations are 
completely neglected and the system is always assumed 
to be described by a single, coherent order parameter. 
Therefore an accurate description of the superfluid to 
Mott insulator quantum phase transition is not possible. 
Nevetheless, qualitative differences exist between mean 
field dynamics in the weak-coupling regime L < 2, where 
the system exhibits full-amplitude Rabi-like oscillations, 
and the strong cou pling ca se L > 2, where self-trapped 
motion is observed |23l M I27L I2I |2<| 

To go beyond the GP approximation, a higher-order 
mean field theory may be formulated by adding to Eq. 
(J2J additional equations of motion for the normal den- 
sity operators djd k , and the anomalous density operators 
dj&k, 




FIG. 2: (color online) Evolution of atomic site populations 
in a two-site system, starting with all population in one site, 
for N = 20, 40, 80 and fixed L = 2. Exact numerical results 
(solid) are compared with the HFB (dotted) and BBR (red 
dashed lines) approximations, as well as to the two approx- 
imations based on the 2PI effective action formalism: 2nd 
order (x's) and l/Af (circles), taken from Fig. 5 in Ref. |23[. 



l -a j a k 



i Jt a jak 



J (a^aj^i + a,f.aj+i + djd k _i + djd k+ i) 



+U(a] 



akajdj + aldjdkCLk 



"2 \ a J a 3 + a kak)(>jk , 



(4) 



J (aja fc _i + a]a fe+ i - a]_ 1 a fe - aj +1 a fc ) 
+U (d^dld k d k - d^d^djdkj ■ (5) 



Taking the expectation values of Eq. J2J and Eqs. J2J- 
and using the HFB Gaussian ansatz, we truncate 
third- and fourth-order moments as: 



(ABC) w (A)(BC) + (B)(AC) 

+ (C)(AB)-2(A)(B)(C) , 
(ABCD) « (AB)(CD) + (AC)(BD) 

+ (AD)(BC)-2(A)(B)(C)(D) 
to obtain the HFB equations: 



(6) 
(7) 



dr 



dr ° k 



+La*a j a j + L (2a, + "1 A" 7 ) , 



(8) 



-2L(\c 



AO 

- A™ 



AV0 



A£J A" 



+ -(a 2 fc + A^) 



2 



a;;,) (2a$ • 



= A 



2L[(|a fc 



A 



_ A" — A" 
1 A j-l,fc A j 



+l,k) 



(9) 
(10) 



(K| 2 + A^)]A 



(4 + A£ fc ) A", 



jk 



*3\ 1 — jj; 
^ + A^) 



jk 



a;, 



for the mean field 



correlation functions A™ fe 



= (dj)/vN and the two-point 



\{a)d k ) 



a*a k ]/N, A% 



[(djd k ) — ajCLk] /N, constituting the reduced single par- 
ticle density matrix. 

We note that the discrete HFB equations JSJ-JTDJl are 
not UV divergent due to the natural momentum cutoff 
imposed by the lattice. However, due to the existence of 
a noncondensate anomalous density, U(l) symmetry is 
broken, in contrast to the gauge-invariant original field 
equations (J2J). U(l) symmetry may be restored for ex- 
ample, by ommitting all anomalous quantities, to obtain 
the Hartree-Fock equations 



■ d A 

1— A 



L(\c 



2A« ) a, 



(11) 



-A% = (A^+A^-A^-A*^) (12) 



-2L[(\a k \ 



AL) ~ (l a j| 2 + A™-)] A" fc . 



Equations (|I1 \ and (|I2|I . conserve separately the conden- 
sate population Wjl 2 an d the noncondensed fraction 
■ AJj. Thus, the HF approximation can not be used to 
account for condensate depletion during the evolution. If 
only the noncondensate anomalous terms are neglected, 
one obtains the HFB-Popov [24j approximation, which 
allows for growth of fluctuations, but conserves the con- 
densate population, so that the total number is not a 
constant of motion. In the following section we construct 
a C/(l) invariant mean-field theory which conserves the 
total number of particles, yet includes dynamical deple- 
tion. 



4 



N=20 



N=40 



N=80 



*0.9 




FIG. 3: (color online) Evolution of the leading eigenvalue 
(above) and single-particle entropy (below) for a two-site sys- 
tem with N — 20, 40, 80 and L — 2. Exact many-body numer- 
ics (solid blue line) is compared with the HFB approximation 
(green dotted line) and the BBR approximation (red dashed 
line) . 



III. THE BOGOLIUBOV BACKREACTION 
EQUATIONS 



Instead of the conventional mean- field approaches, 
based on the site field operators a,j , we construct a mean 
field formalism using the reduced single-particle density 
operator ajafc, treating it as the fundamental quantity. 
We have previously applied this approach to the case of 
a two-site model [26|,|27]]. Here we extend it to the general 
M site case. It is convenient to rewrite the Hamiltonian 
||TJ in terms of the M 2 — 1 traceless operators which gen- 
erate SU(M): 



u h k = aj&k + a\a,j, 1 < k < j < M 

ijM = -i [a]a k - a\a,j J , 1 < k < j < M 



(13) 



Tsl^nj- lhi+i , 1 < I < M - 1 



Since it is easily verified that 

M-l 



1 M-l M 

1 y & + ±n 2 = y 

2 f-f J M f-f 



(14) 



3=1 



3=1 



where n = ^ s ^ ne total particle number, equation 

|JT|| can be rewritten, eliminating c-number terms, as: 



M-l 



h = j J2 fi j+u + j E 

3=1 3=1 



M-l 



(15) 



Using the SU(M) generators we construct a pseu- 



dospin vector operator, 

S = (U21,U 32 , ■ ■ ■ ,«31,&42, ■ • • ,#21, #32, 

. . . , U31,«42, ■ ■ ■ ,Wi,W2, ■ ■ ■ ,WM-l) , (16) 

so that the Hamiltonian (|15|l takes the form: 



.W i M 2 -l 



^E^ + j E § . 

3=1 



2 

4 ^ ~3' 

j=M 2 -M 



(17) 



The Heisenberg equations of motion for the operators Si 
and their products S^S 1 ; then read: 



.d a 



M-l M-l 

3 = 1 fe=l 

M 2 -l M 2 -l 



(18) 



U 
1 



y y <*M +sjs k ), 

j=M 2 -M+l fe=l 
M-l M 2 -l 



dt 



= jy y (c^SkS+c^sA) 



(19) 



i=i fe=i 

M 2 -l M 2 -l 



j=M 2 -M+l fc=l 
M 2 -l A/ 2 -l 



C7 



■j E E c tjSi(S k Sj + S'jS'fe), 



j=M 2 -Af+l fc=l 

where the coefficients Cy are the structure constants of 
the SU(M) group. We note that for M = 2 the Hamilto- 
nian (fTS|> and the dynamical equations l(THl) - (tT^l) reduce 
to the familiar Bloch forms used in Rcfs. |26|,|27(. The 
M -site system is a direct extension of the two- mode case, 
in that hopping terms induce linear Rabi-like oscillations 
in the vw subspace, whereas on-site interactions lead to 
nonlinear phase precession in the uv subspace. 

The reduced single-particle density matrix is obtained 
from the expectation value of S, according to: 



N 



M-l 



\ E 



(20) 



3=1 



where X is a unit matrix of order M and Uj are the 
M x M irreducible representations of the SU(M) gen- 
erators (e.g. Pauli matrices for M — 2, Schwinger ma- 
trices for M — 3 etc.). We will therefore focus on the 
dynamics of the 'hyper-Bloch-vector' S = (S) /27V. The 
lowest-order mean-field approximation replaces the vec- 
tor of operators S by the vector of their expectation val- 
ues S, thus truncating (SiSj) » (Si)(Sj). This results 
in the nonlinear pseudospin-precession form of the GP 
equations, 



S = B(S)®S 



(21) 
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FIG. 4: (color online) Site-populations in a three-site system 
as a function of rescaled time r for N = 20, 40, 80 and fixed 
L = 2. Blue, green, and red lines correspond to 1st, 2nd, 
and 3rd site populations, respectively. Solid lines depict the 
full many-body dynamics, whereas dotted and dashed lines 
correspond to the HFB approximation and the BBR approx- 
imation, respectively. 
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5 10 5 10 5 



T 

FIG. 5: (color online) Leading eigenvalue of the reduced 
single-particle density matrix p and single-particle entropy 
— Tr(plnp), as a function of rescaled time r in a three-site 
system with N = 20, 40, 80 and L — 2. Exact results (solid 
blue lines) are compared to HFB calculations (dotted green 
lines) and BBR calculations (dshed red lines). 



tation values. Taking the expectation values of Eq. 1|18|) 
and Eq. H19[) and truncating 

+(S k )(S i S j )-2(S i )(S j )(S k ) , (23) 
we obtain the BBR equations of motion: 



where 

B(S) = (Bi,5 2 ,...,B M 2-i) , (22) 

with 

r i j = i,..,i(i-i)/2 

Bj ■ = I j = M(M - l)/2 + 1, . . . , M(M - 1) 
[ LSj j = M(M — 1) + 1, . . . , M 2 - 1 

It is readily verified that Eq. I|21l) is exactly equiv- 
alent to the discrete GP equation J2J. In addition to 
the conservation of the total number Tr(p) there exist, 
within GP theory, M — 1 independent constants of the 
motion Tr(p ,n ) with m = 2, . . . , M — 1. For example, 
for M = 2 the GP mean-field theory also conserves the 
single-particle purity Tr(p 2 ), which is just the length of 
the three-dimensional Bloch vector. Deviations from this 
classical field theory, due to interparticle entanglement 
and loss of single particle coherence, will show up as a 
reduction in these classically conserved quantities. 

The BBR approximation is obtained by going one level 
deeper in the hierarchy of dynamical equations for expec- 



, M-1M 2 -1 

j=\ k=l 

M 2 -l M 2 -l 

+ L E E 4(SA + A, fe ) (24) 

j=M 2 -M+l k=l 
, A/-1M 2 -1 

V Aj ' = E E (c&Att + cfcAtt) (25) 
T j=i k=i 

M 2 -l M 2 -l 

+ L E E 4(A l3 S k + A lk S 3 ) 

j=M 2 -M+l k=l 
M 2 -l M 2 -l 
+i E E efj (&ijSk + &ikSj) , 

j=M 2 -M+l k=l 

where Sj = %)- and A jk = <5j5fc+5 4 fc ^ ) ~ 2Sj5fc . In the 
following section we compare the accuracy of the BBR 
approximation with respect to GP, HFB, and 2PI effec- 
tive action. 
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FIG. 6: (color online) Site-populations in a four-site system 
as a function of rescaled time r for N = 20, 40, 80 and fixed 
L = 2. Blue, green, red, and cyan lines correspond to 1st, 
2nd, 3rd, and 4th site populations, respectively. Solid lines de- 
pict the full many-body dynamics, whereas dotted and dashed 
lines correspond to the HFB approximation and the BBR ap- 
proximation, respectively. 



IV. NUMERICAL RESULTS 

In order to test the accuracy of the BBR approxima- 
tion compared to other methods, we carried out exact 
numerical calculations for limited numbers of particles 
and sites (up to N — 100 particles and M — 5 lattice 
sites). The Hamiltonian was represented in terms of 
site-number states and the iV-body Schrodinger equation 
was solved numerically, as in Refs. |U Efj l27| . These 
many-body results were then compared with BBR mean- 
field calculations, as well as with GP, HFB and variants 
of the 2PI effective action method. 

In Fig Q] the evolution of fractional population differ- 
ence for a hundred particles in two-sites, is plotted for 
various values of the coupling parameter L. Within the 
GP mean-field theory, full-amplitude Rabi-like oscilla- 
tions are predicted in the linear regime with L < 2 (Fig 
Ha)). As the transition is made to the strong-coupling 
regime, the oscillation becomes increasingly more non- 
linear, until when L > 4 macroscopic self-trapping is 
attained (Figs. E^b)-^d)). The value of L = 4 is par- 



ticularly interesting because for this coupling a trajec- 
tory starting from a single-populated site becomes dy- 
namically unstable when site-populations equilibrate. In 
previous work we have shown that this dynamical in- 
stability serves as a quantum- noise amplifier |26l |27| , so 
that the growth of the deviation of a quantum trajec- 
tory from the corresponding GP prediction is initially 
exponential, leading to a log(l/N) slow convergence of 
the many-body quantum-field results to the classical GP 
prediction. Thus, while the naive expectation would be 
that quantum fluctuations would simply grow with the 
coupling parameter L, their role is in fact maximized for 
L = 4, as evident in Fig. [Jb). It is clear from Fig. Q] 
that the BBR approximation gives a better description 
of the ensuing quantum dynamics, for longer timescales, 
than HFB does. 

Convergence of various approximations with increas- 
ing number of particles is demonstrated in Fig. [21 where 
the two-sites population dynamics is plotted for increas- 
ing particle numbers, keeping a fixed coupling value of 
L = 2. In addition to the exact, BBR, and HFB results, 
we also plot two calculations based on the 2PI effective 
action approach, taken from Fig. 5 of Ref. |2^| (our exact 
and HFB results exactly coincide with the corresponding 
lines in that figure). Here too, the BBR approximation 
(red dashed lines) gives a more accurate description of 
the dynamics than any of the other methods, attaining 
a nearly perfect convergence in the given time-frame for 
N = 80 particles. In comparison, standard HFB fails to 
depict the damping of coherent oscillations, whereas the 
2PI effective action methods tend to overdamp. We note, 
that in terms of formalistic complexity alone, the BBR 
approximation is far simpler than the noninstantaneous 
integrodifferential equations used in the 2PI effective ac- 
tion methods [2j|. In fact, it is even simpler than HFB, 
in that only normal quantitities are involved, giving a to- 
tal of nine equations for two sites, as opposed to fifteen 
in HFB. 

Dynamical condensate depletion is also well-depicted 
by the BBR approximation. In Fig. [21 we plot the evo- 
lution of the leading eigenvalue of the reduced single- 
particle density matrix p and the single-particle von- 
Neumann entropy — Tr(plnp), corresponding to the pop- 
ulation dynamics of Fig. [21 While HFB calculations seem 
to give an abrupt deviation of the predicted condensate 
fraction from its exact value, the BBR results converge 
well, giving a reasonably accurate description of BEC de- 
pletion. 

The same qualitative behavior carries over to systems 
with more than two sites. In Fig.0]and Fig. [SI population 
dynamics and condensate depletion are shown for a three- 
sites system with N = 20,40,80 particles. Similarly to 
the two-sites case, the BBR approximation constitutes a 
significant improvement over the HFB approach, giving 
a better description of populations as well as coherences. 
The same is also true for the four-sites case shown in Fig. 

El 

The faster convergence of BBR as compared with the 
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log(N) log(N) 

FIG. 7: (color online) Characteristic times at which the Carte- 
sian distance between the exact Bloch vector and its HFB 
(green, x's) and BBR (blue, circles) approximants, reaches a 
predtermined threshold, as a function of N for L = 1 (a), 
L = 4 (b), L = 6 (c), and L = 10 (d). The break-threshold is 
set to 0.2 in (a)-(c), and to 0.05 in (d). 

HFB appoximation is illustrated in Fig. [7| where charac- 
teristic breaktimes of the two approximations in a two- 
sites calculation, are plotted as a function of the total 
number of particles N. As anticipated, breaktimes grow 
as y/~N when the classical dynamics is regular JZ^EHTtO 
and as log A when the classical trajectory hits the dy- 
namical instability Qd) . The BBR calculations give con- 
sistently longer breaktimes, with a more regular conver- 
gence pattern. 

V. DISCUSSION 

The rich regime of strongly correlated many body 
physics, which ultracold atom experiments are now be- 
ginning to probe, will surely not be fully conquered by 
any simple hierarchy truncation scheme such as BBR. 
Nor does BBR offer anything like an exact solution even 
to the problems to which we have applied it in this 
paper; its improvements over its rivals are incremental 
rather than revolutionary. On the other hand it should 
be born in mind that incremental improvements in the- 
ory are more significant in the context of ultracold gases 
than in traditional condensed matter, because in the 
new atomic systems samples are precisely characterized, 
controlled, and measured, and relevant microphysics is 
clearly known. It is perfectly plausible in these systems 
that we may come to learn important qualitative princi- 
ples from experimental discrepancies on the few percent 



level. 

The merits of BBR that we would like to emphasize, 
along with its very reasonable level of accuracy, are its 
simplicity and its direct relation to experimental reality. 
It involves only quantities which are directly observed 
in single- and two-particle number-conserving measure- 
ments, and it respects the fact that in current quantum 
gas laboratories atoms are neither created nor destroyed. 
And it is conceptually and computationally straightfor- 
ward. 

In one sense it is of course conceptually all too straight- 
forward: like all hierarchy truncation schemes since 
Boltzmann's, it is an uncontrolled approximation, whose 
accuracy is therefore arguably as much a puzzle as it is 
a solution. Insofar as truncating a hierarchy at a deeper 
level is grounds for expecting higher accuracy, however, 
the advantage of BBR is clear: it is a truncation at fourth 
order in field operators, compared to only second order 
for HFB. Deeper level truncation often involves prolifer- 
ation of terms, to the point of sharply diminishing re- 
turns in accuracy versus effort; but BBR avoids this, and 
manages to use fewer equations than HFB, because it 
eliminates all anomalous terms. 

And this leads us to conclude by indicating some of the 
potential future applications of the results of this paper. 

Why do hierarchy truncations often work as well as 
they do? What determines the best way of truncating 
a hierarchy? These are questions that have been raised 
ever since Boltzmann's stosszahlansatz produced the ar- 
row of time, but they have yet to be fully answered. 
With current experimental capabilities for precise and 
controlled measurements on ultracold gases, introducing 
a physically motivated alternative truncation scheme, as 
this paper has done, may contribute to new progress on 
these questions. 

Finally, another conceptual merit of BBR is that be- 
cause it is based on the single particle density matrix, 
rather than just the macroscopic wave function, it makes 
such a conceptually important quantity as the single par- 
ticle entropy - the entropy of Boltzmann - a basic in- 
gredient in the theory, rather than a perturbative af- 
terthought. Rethinking entropy, heretofore mainly in the 
context of quantum information and computation theory, 
is another major thrust of current physics; the alterna- 
tive viewpoint offered by BBR may potentially be of some 
value in a broader conception of this project. 
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